treetop: A Shiny-based application and R package for extracting forest information from LiDAR data for ecologists and conservationists

1.1 Background and model foundation

The use of LiDAR for ITD and crown delineation is promising; however, it requires specific training in computational methods for practitioners to derive reliable forest information from LiDAR datasets (Valbuena et al., 2014). A LiDAR-derived CHM can be used to detect individual trees, delineate their crowns and subsequently estimate biophysical attributes such as biomass, stem volume and gap fraction (Silva et al., 2016; Stark et al., 2012). The CHM is defined as the spatially explicit height of the tallest vegetation and is essentially computed by subtracting the digital terrain model (DTM) from the digital surface model (DSM) representing the scanned vegetation (Figure 1). The simplest method to detect individual trees on the LiDAR-derived CHM is by using the local maximum (LM) algorithm. LM methods assume that local height maxima in the CHM represent tree tops (Korpela et al., 2006), and the method is a relatively straightforward technique that utilizes two major parameters: a smoothing parameter, which is also often mentioned as smoothing window size (SWS), and a fixed window size (FWS) for tree detection (Silva et al., 2016). As FWS increases, the number of trees detected decreases (Mohan et al., 2017). Application of smoothing filters helps eliminate spurious local maxima caused by, say, large spreading tree branches, and thereby limits the number of local maxima detected, increasing the accuracy of the algorithm (Lindberg & Hollaus, 2012). Once, individual trees are detected, their crown boundaries can be delimited using the Voronoi tessellation-based algorithm developed by Silva et al. (2016) (Figure 2). The algorithm uses maximum crown factor and exclusion parameters, both ranging from 0 to 1, to define the crown boundaries on the LiDAR-CHM (Silva et al., 2016).

treetop is both an interactive online tool as well as a standalone R package (Silva et al., 2021) that puts complex ITD and crown delineation procedures easily within reach of any non-specialized user, allowing the visualization of forest information promptly and the ability to download generated results in a variety of file formats. The tool was developed using RStudio, a free and open-source integrated development environment for R (R Core Team, 2018), and the Shiny package in R (Chang et al., 2019). The online version of treetop is publicly available on the shinyapps.io, a service platform for hosting Shiny web apps, at https://carlosasilva.shinyapps.io/weblidar-treetop/ and has limited size for input (30 Megabytes). The R package version of treetop can run natively on a user's computer and can be downloaded via the Comprehensive R Archive Network (CRAN) (https://cran.r-project.org/web/packages/treetop/index.html). In addition, it can be installed directly in R install.packages (treetop) and launched by the single command line launchApp(). The R package launches a native Shiny-based app identical to the treetop web-based tool.

We used the LM, SWS and FWS algorithms for ITD (Figure 2) along with three CHM filter options and a Voronoi tessellation-based algorithm for crown delineation. This algorithm and ITD accuracy have been already validated in Silva et al. (2016) and Leite et al. (2020). Optionally, the user can estimate tree crown width using a list of equations available in the application (Popescu et al., 2014) or even by a user-defined custom equation.

1.2 System design and features

The requested data input for ITD through the treetop application is a LiDAR-derived CHM in raster format, and thus CHM-derived from other remote sensing sources can also be used. The R package version of treetop does not have any data size limitations, thus allowing users to process larger CHM rasters. The treetop tool is capable of providing crown change estimates from multitemporal data as highlighted in our first case study. The treetop application graphical user interface (GUI) has five major panels: (a) the settings menu; (b) the summary of the LiDAR metrics; (c) the CHM profile (frequency distribution and box and whiskers plots) or spatial pattern analysis; (d) the individually detected trees represented over the CHM or Lorenz curve analysis of tree size heterogeneity; and (e) an interactive 3D viewer of the CHM or a reconstruction of individually detected trees (Table 1; Figure 3). Furthermore, there are eight buttons to download the generated results, depending on the users' choices, which include: (a) metrics summarizing the distribution of LiDAR heights in the CHM; (b) a figure of the CHM profile, which includes the frequency distribution and box plots; (c) a figure showing the detected trees and crowns over the CHM; (d) the LiDAR-derived trees with the tree location and crown width estimates; (e) two buttons for downloading the resulting tree positions and delineated crowns in shapefile format; (f) spatial pattern analysis of detected trees; and (g) Lorenz curve analysis of differences in height among the detected trees.

2.1 Forest disturbance detection using airborne and UAV-LiDAR multitemporal data and the treetop application

Hurricanes cause extensive damage to forest resources, in addition to the loss of life and other properties. The ecological effects of hurricanes on forest structure and function are not well-understood. They can be complex, including both harmful and beneficial effects. Assessing damages and monitoring forest structural changes following hurricanes is the first step to understanding complex ecological effects of hurricanes on forest and develop forest management approaches aimed at increasing forest resistance and resilience (Shiels et al., ).

Ecologically, hurricanes are part of a suite of natural disturbances that have historically shaped structure and function of forest ecosystems in the Gulf and Atlantic Coastal Plains of the United States. On 10 October 2018, a powerful hurricane named Michael made landfall in the Florida Panhandle (Klotzbach, 2018). The present case study was conducted in the Apalachicola National Forest in Florida Panhandle following hurricane Michael, which severely damaged large parts of this forest. We made a quick assessment of damage and disturbance on sites impacted by Hurricane Michael using LiDAR data and the treetop tool. Specifically, we used LiDAR data that were collected immediately prior to Hurricane Michael from an aircraft-borne LiDAR system. A UAV-LiDAR system (GatorEye system, see Supporting Information) collected data at the same site 2 months after the hurricane which was 6 months after the initial airborne flight. CHMs (0.5 m resolution) were computed using LAStools (Isenburg, 2018) and, for the purposes of this case study, we selected a 50 × 50 m sample plot in which individual trees were detected and delineated on the CHMs using FWS = 3 × 3 and SWS = 3 × 3, and MaxCrown and Exclusion parameters were set to 0.6 and 0.3 respectively (Figures S1 and S2).

Hurricane impacts were visually apparent between the collection dates in our sample site. The 32 trees initially detected were reduced to 24 trees following the hurricane. We found 38.5% reduction in canopy coverage with large increases in gaps and crown reductions in some remaining standing trees (Figure 4). No decrease in maximum tree height was detected within the plot; however, mean tree height decreased from 26.8 to 23 m, indicating that mortality of the canopy trees possibly occurred. Height distributions (Figure 5) can be used to calculate vertical leaf area density provided sufficient calibration data (Stark et al., 2012). Even without calibration, the assessment of these vertical profiles reveals a large increase in LiDAR echoes below 8 ms above the ground which is consistent with an increase in downed woody debris following the hurricane. However, for this small sample plot, a two-sample Kolmogorov–Smirnov (KS) test (Wayne, 1990) showed that the tree height distribution pre- and post-hurricane disturbance (Figure 5) did not significantly differ (KS = 0.23, p-value = 0.47 and α = 0.05).

2.2 Assessment of structural heterogeneity in forest plantations

The treetop tool was used to access forest structural heterogeneity in a 3-year-old fast-growing Eucalyptus grandis forest plantation (e.g. Silva et al., 2014) and in an 8-year-old seminal (seed) plantation of Pinus taeda (Figure 6) in southern Brazil (Klabin SA company, Paraná State, Brazil; Silva et al., 2017). LiDAR data were collected with pulse density of 4 pulses/m² and CHM (0.5 m of resolution) was created using LAStools (Isenburg, 2018). Individual trees were detected and delineated in both plantation sites using treetop tool in a plot area of 0.25 ha (50 × 50 m²). ITD and crown delineation were performed using FWS = 3 × 3 and SWS = 3 × 3, and MaxCrown and Exclusion parameters were set to 0.6 and 0.3 respectively (Figures S3 and S4).

The structural heterogeneity of forests is typically assessed by a measurement of variability and a measurement of asymmetry (Knox et al., 1989; Valbuena et al., 2017). In our treetop tool, we chose to evaluate variability and asymmetry in the height of detected trees, because height above the ground is the variable directly measured by the LiDAR. Except for MacArthur and MacArthur's (1961) foliage height diversity (FHD), structural indices are more typically based on tree bole diameters because in the field they are more directly measurable than heights. One such measure of heterogeneity based on tree heights is the PH³50 index (Hentz et al., 2018). To clarify the role of PH³50 in evaluating forest heterogeneity, we computed Lorenz curves from the cubic power of tree heights, which has been observed to relate to Lorenz curves obtained from tree diameters (Hakamada et al., 2015). Lorenz plots for the two forest areas—(a) the E. grandis and (b) P. taeda plantations (Figure 7)—both showed a high degree of homogeneity in their structure, which can be observed by the low amplitude of their Lorenz curves, close to the diagonal (absolute equality). The Gini coefficient obtained from both of stands, GC = 0.17 and 0.11, respectively, are reflective of the Lorenz curve results (Weiner & Thomas, 1986). The upper part of the shaded area in Figure 7 shows the position of a Lorenz curve with a maximum value for MacArthur and MacArthur's (1961) FHD, as demonstrated by Valbuena et al. (2012) depiction of maximum entropy over the Lorenz plot.

More apparent differences between plantations can be observed when considering the measures of asymmetry. Damgaard and Weiner (2000) envisaged a method to evaluate the asymmetry of a Lorenz curve by comparing the position of its inflection point (a crossed circle in Figure 7) against the axis of symmetry of a Lorenz plot (dashed line in Figure 7). Complete symmetry is denoted by an index value of S = 1 (Damgaard & Weiner, 2000). Our case study plots showed values of 1.03 and 1.18, respectively, which reveals the P. taeda forests to have a more heterogeneous structure. In these two case studies, we observed values of PH³ = 0.55 and 0.64, and PH³50 = 0.58 and 0.61 respectively (in Figure 7, PH³ is denoted by a crossed circle and PH³50 by a red arrow). These are just a few examples of how ecologists, forest science experts, foresters and other practitioners can benefit from the outputs automatically yielded by the Web-LiDAR treetop application and the R treetop package.

2.3 Spatial patterns of individual trees in tropical forest

In our third case study, the treetop tool was used for ITD and crown delineation within a plot of 2.25 ha in a tropical forest site in the Western Brazilian Amazon to analyse the spatial distribution of trees by using Ripley's K functions (Ripley, 1976). The LiDAR data were collected at the Adolfo Ducke Forest Reserve in Manaus, state of Amazonas, Brazil (dos-Santos & Keller, 2015) with a pulse density of 22.7 points/m² in 2012. CHMs (0.5 × 0.5 m) were created using LAStools (Isenburg, 2018) for ITD and crown boundary delineations (Silva et al., 2016). ITD was performed using FWS = 7 × 7 and SWS = 7 × 7, and MaxCrown and Exclusion parameters were set to 0.4 and 0.7 respectively (Figures S5). The treetop detected 75 trees with mean (± SD) height of 33.47 m (± 4.93) and a mean crown area of 213.79 m² (± 86.96) (Figure 8).

Ripley's K and L functions (Besag, 1977; Ripley, 1976) and the Clark–Evans index (Clark & Evans, 1954) were used to examine spatial patterns followed by the detected trees in the 2.25-ha plot in the Adolfo Ducke Forest Reserve. Ripley's function calculates distances (r) between individual pairs of trees and compares them to what would be expected from a purely random (Poisson) point process. The empirical K function observed from spatial distribution of trees (Kobs (r); black solid lines in Figure 9) is compared against the theoretical K function (Ktheo (r) = πr²; red dashed line in Figure 9). The grey-shaded areas in Figure 9 are confidence intervals for the random process developed by Monte Carlo simulations (Baddeley et al., ). Random variability and large ranges of values in the point intensity (y-axis) makes it difficult to analyse the differences between both curves (Figure 9a). For this reason, Besag (1977) proposed the L function, as a transformation of the K function, Lobs (r) = 𝐿(𝑟)/𝜋⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯√$\sqrt{L\left.r\right)/\pi }$, where the theoretical L function becomes Ltheo (r) = r (Figure 9b). Finally, a further transformation into the 𝐿̂$\widehat{L}$ function, obtained after subtracting the distance r from the L values, which makes Ripley's function and the confidence intervals for a random process easiest to visualize and analyse (Baddeley et al., ). Clark and Evans (1954) index (R) summarizes the information of Ripley's function into a single indicator with values R = 1 representing a random process, R < 1 a mechanism of attraction resulting in aggregated spatial patterns and R > 1 a mechanism of inhibition which results in spatial regularity (Figure 9c).

The results obtained in the Adolfo Ducke Forest Reserve plot with the treetop application showed signs of mechanisms inhibiting the growth of large trees close to one another (Figure 9). This is revealed by the empirically observed Ripley's K function laying below the confidence intervals for random process for distances around r = 4–15 m. The overall value of Clark-Evans index R = 1.48; from an ecological viewpoint, this indicates the mechanism of size-asymmetric competition between these trees which inhibit the growth of a tree in proximity of another (Weiner, 1990), which is in agreement with the natural growth conditions of this forest (DeMalach et al., 2016; Freckleton & Watkinson, 2001).